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ABSTRACT 

A remarkable phenomenon in turbulent flows is the 
spontaneous emergence of coherent large spatial scale 
zonal jets. Geophysical examples of this phenomenon 
include the Jovian banded winds and the Earth's po- 
lar front jet. In this work a comprehensive theory 
for the interaction of jets with turbulence, Stochastic 
Structural Stability Theory, is applied to the problem 
of understanding the formation and maintenance of 
the zonal jets that are crucial for enhancing plasma 
confinement in fusion devices. 

1. Introduction 

Coherent jets that are not forced at the 
jet scale are often observed in turbulent flows 
with a familiar geophysical example b eing the 
zonal winds of the gaseous planets (jlngersoU 
199Ci) . This phenomenon of spontaneous jet for- 
mation in turbulence has been extensively inves- 
tigated in observational and theoretical studie s 
( Balk et al.lll990l: |Panettalll993t [ValHs and Maltru 



nl 



1^9931; 'Cho and Polvanil Il996' 'Huang and Robinsoni 
"998- Farrell and loannou 2003 2007: Diamond et al 
2005t IConnaughton et al. 20091) as well as in labo- 



Fuiisawa et al. 


2008: Itoh et al. 2007a,b: Read et al] 


20071 Mazzucato et al.l Il996l iHolland et all l2006l) 



The mechanism by which these zonal flows 
form and are maintained is systematic organi- 
zation of upgradient eddy momentum flux in 
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which the transfer of momentum occurs directly 
from the eddy field to the zonal fiow without 
passing through intermediate scales, in contrast 
to the prediction of theories based on two di 



mens io nal turbulence cascad e s (iNozawa and Yoden 
| l997 : Huang and Robinson 1998t Ingersol et al 



'2004': 'Salvk et al."2006l: iKitamura and Ishiokal 12007 
Diamond et al. 2005). 

Excitation of the eddies that give rise to zonal jets 
in turbulence can be traced either to predominantly 
external processes such as convection, as in the case 
of the Jovian jets, or to predominantly internal pro- 
cesses such as baroclinic growth, as in the Earth's po- 
lar front jet. However, maintenance of turbulence in a 
given flow is usually due to a combination of external 
and internal processes as for instance latent heat re- 
lease associated with cumulus clouds injects external 
potential vorticity perturbations into the baroclinic 
turbulence of the polar front jet. 

Because of their self-regulating nature and inter- 
dependence drift wave turbulence and zonal flows 
behave as a single drift wave - zonal flo w system 
(hereafter DW-ZF) (jDiamond et al.ll2005f ). In this 



system the drift wave perturbations arise from the 
internal instability of the imposed density gradi- 
ent, from sources external to the intrinsic dynam- 
ics of the drift waves and at a given scale from 
transfer between scales by the internal quadratic 
nonlinear advection. Because these processes pro- 
duce perturbations with short time and space scales 
compared to the time and space scale of the jet, 
the associated eddy dynamics can be simulated us- 
ing a Stochastic Turbulence Model (STM) in which 
the nonlinear scattering and extrinsic excitation are 
modeled as stochastic ( Farrell and loannoul 1993a . 
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: Newman et al. 1997; 


Zhane and Held 19991 DelSold 


2004). The STM pro- 



vides an analytic method to obtain the dynamics of 
the quadratic statistics of a turbulent eddy field as- 
sociated with a given jet structure. Coupling a time 
dependent STM to an evolution equation for the jet 
produces a dynamical system for the co-evolution 
of the jet and the self-consistent quadratic statis- 
tics of its associated turbulence; this is the method 
of Stochastic Structural Stability Theory (hereafter 
SSST). The SSST system can be interpreted as the 
dynamics of the ensemble mean jet and the ensemble 
mean associated turbulence in which the turbulence 
is modeled by the ensemble mean perturbation lin- 
ear dynamics with a stochastic approximation for the 
non- linear dynamics. The solution for the eddy field 
is in terms of a covariance matrix from which can 
be obtained the Gaussian probability density func- 
tion approximation for the variance and quadratic 
fluxes of the turbulence. The solution trajectory of 
the SSST equations often converge to a fixed point 
state of balance between the turbulence and the jet; 
how ever, limit cycles and chaotic solutions also oc- 
cur (jFarrell and loannoullioosl . l2007i I200I l2009al lbh. 



Chaotic trajectories of the SSST system correspond 
not to chaos of an individual turbulent state tra- 
jectory, which typically would be associated with a 
fixed point of the SSST system, but rather to chaos 
of the ensemble mean turbulent state itself. A fa- 
miliar example of this type of chaos is the irregu- 
lar bursting behavior s een in drift wave turbulence 



([Mazzucato et al.lll996[ ). Conceptually it is useful to 
view the SSST system as a computationally tractable 
approximation to a deterministically initialized Liou- 
ville system for the associated flow. 

Interaction between the zonal flow and its consis- 
tent field of turbulent eddies is nonlinear and can 
support multiple equilibrium states. In many cases 
these equilibrium states arise from and can be traced 
by continuation in a parameter to a bifurcation of 
the coupled DW-ZF system. In the case of the equi- 
librium state with no mean density gradient and no 
initial zonal flow the zonal flow forming bifurcation 
arises as a function of a parameter controlling tur- 
bulence intensity as an emergent instability of the 
SSST system intrinsic to the interaction between the 
zonal flow and the turbulence. One may think of a 
perturbation zonal flow organizing the surrounding 
turbulence to produce a momentum flux divergence 
that amplifies that perturbation zonal flow. The par- 
ticular perturbation zonal flow structure that orga- 
nizes the turbulence to exactly amplify its own struc- 
ture is obtained as an eigenfunction of the perturba- 
tion SSST system linearized about a marginally sta- 



ble SSST equilibrium. This instability equilibrates at 
finite amplitude and this finite amplitude SSST equi- 
librium, consisting of the zonal flow and associated 
consistent eddy field, can be connected by continu- 
ation in an appropriate parameter, such as the den- 
sity gradient, to nearby finite amplitude equilibrium 
states. 

In addition to simply continued equilibria there 
also exist equilibria that are isolated to variation 
of a given parameter as for instance a strong zonal 
flow equilibrium exists at a moderate density gradi- 
ent and turbulence intensity that can not be con- 
nected by continuation starting from a weak zonal 
flow equilibrium at a small density gradient. How- 
ever, external turbulence excitation can be used as a 
control parameter to promote the system to such an 
isolated equilibrium state. In addition to parameter 
control we may also perturb the zonal flow to promote 
the system to an isolated equilibrium state. Pro- 
moting the DW-ZF system to different regime states 
by parameter control is analogous to instigating a 
laminar/turbulent transition in shear flow turbulence 
where the Reynolds number is the control parameter. 

An equilibrium state of balance between a zonal 
flow and its associated field of drift wave turbulence 
requires that the momentum flux divergence arising 
from the turbulence precisely balance the zonal flow 
momentum loss to friction, if any. The requirement 
of a precise balance between zonal flow forcing and 
dissipation, if any, is far more demanding than that 
the shear associated with the jet simply suppress the 
turbulence while the turbulence during the suppres- 
sion process produces up-gradient momentum flux. 
The remarkable fact is that the turbulence, which 
depends on the zonal flow, and the zonal flow, which 
depends on the turbulence, mutually adjust to pro- 
duce balanced states. Having the SSST analytic dy- 
namics of the DW-ZF system allows us to predict pa- 
rameter values for which robust equilibrium DW-ZF 
regimes are maintained, to predict parameters values 
for which time dependent periodic and chaotic DW- 
ZF regimes occur, to predict transition between these 
regimes when two regimes exist at the same param- 
eter values, and ultimately to predict the breakdown 
of the zonal flow regime. 

Closer inspection of the density transport mecha- 
nism reveals that the observed and simulated DW-ZF 
equilibrium jet density transport suppression can not 
be u nderstood using the concept of effective diffu- 



sion ( Sanchez et al. In effective diffusion the- 



ory it is assumed that transport of a passive scalar 
is proportional to the scalar gradient with coefficient 
Deff — vl in which v and I are the characteristic 
velocity and spatial correlation scales of the turbu- 
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lence. Transport can vary either due to changes in 
the characteristic velocity or in the eddy correlation 
scale. In this work we solve for the correlation be- 
tween velocity and density fluctuations directly re- 
vealing turbulent transport both up and down the 
mean gradie nt, in agreement with observations an d 
simulations (|Shats et al.ll2000l : [Holland et al.ll2006f) . 
and implying the density transport process in drift 
wave turbulence is not diffusive in nature. Instead we 
find that large scale coherent structures rather than 
small scale eddy diffusion are responsible for density 
transport IjBos et al.ll2008r ). 

Closer inspection of the dynamics of the interac- 
tion between perturbations and zonal flows reveals 
that understanding reduction of turbulence variance 
by zonal flows through the concept of shear suppres- 
sion by zonal flow advection is incomplete. Shear 
suppression has roots in WKB theory and the con- 
cept of a continuous spectrum of advected harmonic 
waves. However, to properly understand perturba- 
tion dynamics in jets a full wave solution must be ob- 
tained because the perturbation dynamics supports a 
complete set of large scale coherent modes that are 
in general not orthogonal and among which exists a 
subset that is potentially unstable. Interaction be- 
tween the zonal jet and the eddy field systematically 
stabilizes these modes ( loannou and LindzenI Il986t 
James! 1987 : Roe and Lindzen Il996h during the es- 
tablishment of a statistical equilibrium. Moreover, 
the non-normal equilibrium jet dynamics supports 
a subset of stable structures that produce robust 
growth under internally and externally imposed exci- 
tation. These Stochastic O ptimal (SO) perturbations 
( Farrell and Ioannoulll996l ) comprise a small subset 
of structures but these are the structures responsible 
for growth of perturbations due to interaction with 
the zonal shear and density gradient. Using SSST 
we solve for the complete normal mode eigenstruc- 
ture of the equilibrium jet as well as the SO and 
EOF (Karhunen-Loeve) decomposition of the ensem- 
ble mean turbulence variance and cross variance in 
the velocity and density fields. This analysis pro- 
vides full information on the structure and dynam- 
ics of the perturbations responsible for producing the 
turbulence variance and fluxes. 

The mechanism of jet formation in plasmas can be 
studied for turbulence arising from external, internal, 
or a combination of sources. The Charney-Hasegawa- 
Mima (C-H-M) equation provides the simplest model 
system as it uses only external turbulence excita- 
tion. Zonal jet formation in this model is identical 
to t hat in the equivalent baro tropic vorticity equa- 
tion (^rreir^iJoanno3[2QQ3)- However, because in 
the DW-ZF problem there exists an internal instabil- 



ity associated with the density gradient this problem 
is more comprehensively modeled using the modified 
Hasagawa-Wakatani (H-W) equations which describe 
plasma dynamics in a 2D slab model. These equa- 
tions are similar, altho ugh not identical, to the b aro- 
clinic two-layer model (jFarrell and loannoullioOSh . In 
this work we use the H-W equations to study DW-ZF 
dynamics. 



The SSST equations incorporate a stochastic tur- 
bulence model but these equations are themselves de- 
terministic and autonomous with dependent variables 
the zonal flow and the ensemble mean covariance of 
the turbulence. It follows that the perspective on 
stability provided by these equations differs from the 
more familiar perspective based on the perturbation 
stability of the zonal flow. In fact the primary bifur- 
cation in these equations has no counterpart in zonal 
flow stability analysis; it is rather a cooperative in- 
stability in which the perturbation zonal flow orga- 
nizes the background turbulence to produce flux di- 
vergences configured to amplify the jet leading to an 
emergent turbulence-zonal flow instability that need 
not coincide with perturbation instability of the jet. 
The SSST equations are the nonlinear ensemble mean 
dynamics for the DW-ZF flow system and this system 
in many cases supports equilibration of the emergent 
jets and their consistent turbulence fields at finite am- 
plitude. These finite amplitude equilibria in turn lose 
structural stability as a function of parameters and 
this instability is associated either with bifurcation to 
another equilibrium or to loss of a stable equilibrium 
state. While it is true that loss of modal jet stability 
by an equilibrium state as a function of a parameter 
also implies loss of structural stability, the converse 
is not true. For this reason bounds on zonal jet am- 
plitude based on modal instability of the jet are not 
tight and can often be improved by analysis of the 
structural stability of the jet. 



A gradient driven flow with constant density gra- 
dient is assumed in the examples below for simplicity 
although the particle flux is calculated and could be 
used with an appropriate density gradient forcing pa- 
rameterization to obtain equilibria in which the den- 
sity gradient participates in the equilibration. How- 
ever, as equilibrium is approached the fluxes are typi- 
cally suppressed implying long time scales for changes 
in the equilibrium density gradient by flux divergence 
and the likelihood that external driving mechanisms 
dominate density gradient variation. 
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2 . Formulat ion 

a. The Hasegawa- Wakatani drift wave turbulence 
equations 

We use the modified H asegawa- Wakatani (H-W) 
equations (jNumata et al ] [20071) . These equations 
model the turbulence of the edge region of a toka- 
mak plasma with fractional density decreasing in the 
radial direction, x, at a constant rate k, so that 
n{x) = uqc^'^^, and in a constant background mag- 
netic field B — BqZ in the toroidal, z, direction. The 
H-W equations govern the dynamics of the electro- 
static potential ecji/Tf, and the ion density n/no in a 
cartesian approximation of the radial-poloidal, x — y, 
plane. 

The ion vorticity, Q = A(^, and the density fluctu- 
ations, n, solve following iNumata et al.l (j2007l ): 



dtC + J[(^X]=a{(t>'-n')^vA(^, (la) 
dtn + J[cj), n] = a{(j)' - n) - ndycj) + vAQ (lb) 

with Jacobian J[/, g] = {d^f){dyg)-{dyf){d^g). The 
fields are decomposed into zonal means and depar- 
tures from zonal means: 



, n = n + n', 



(2) 



with the zonal mean, denoted by a bar, defined as 
the mean in the poloidal direction y: 



(3) 



/ = L^M f{x,y,t)dy 



The flow velocities are: 



-dy4> , V ^ dx 



(4) 



The parameter a controls the strength of the elec- 
tron resistivity that couples the electrostatic field 
with the ion density perturbations. For a = equa- 
tion pap for the electrostatic potential corresponds 
to the hydrodynamic 2D vorticity equation while the 
density equation (|lb[) corresponds to the advection- 
diffusion of n' as a passive scalar in the presence of 
a mean fractional radial density gradient —n. In the 
limit a — > 00 the density and electrostatic field cou- 
ple rigidl y and obey the Charney-H asegawa-Mima 
equation (jHasegawa and Mimal 119781 ). The dynam- 
ics of this equation, which governs the formation of 
zonal fiows in both the GFD and the plasma context, 
has been studie d in recent theoretical work on zonal 
flow generation (Balk et al.''l990':'Connaught on et al.l 
[2009; FarreU and loannou 2003, 2007). Hereafter, 
we treat the more general quasi-adiabatic case with 
a — 1, and allow for instability by including an ion 



density gradient, k, which will be treated as a vari- 
able parameter. 

In the nondimensionalization of the equations 
lengths are scaled by the Larmor radius ps — 
\jTf,/miU!~j^ and time by the electron cyclotron fre- 
quency ujci = eBo/rrii. A typical Larmor radius, 
Ps = I mm, is obtained for a magnetic field 1 T and 
electron temperature Tg = 95.6 eV; also for these val- 
ues LLi~^ = 10~^ s/rad and the corresponding velocity 
scale Psi^ci is 95.6 km/ s. The channel is taken doubly 
periodic in both x and y. 

The zonal average of (fTa|) gives the equation for the 
zonal jet: 

dtv = -u'C - rraV . (5) 

Where v = D(f> and D = dx- The zonal fiow is 
damped linearly at the mean coUisional damping rate, 
rm, which will typically be taken to be r„i = 10^^ al- 
though we will also present results in the coUisionless 
limit, rm = 0. 

The nonzonal components obey the equations: 



dtC = -vdyC + {D^v)dy(t>' -f a 
dtn' = —vdyu' — Kdy4>' + a{4. 

with nonlinear scattering term: 

F{f) = -dx{u' r -Vf) - dy{v' f 



- n') + z/AC' + F(C') 
(6a) 

n') + i^An' + F{n') . 
(6b) 



v'f) ■ (7) 



These equations can sustain turbulence without 
external forcing due to the radial density flux, 
u'n', in the presenc e of the mean density gradient 
(|Numata et al.ll2007l ). We now briefly review the en- 
ergetics of these equations. The total energy, E, is 
the sum of the zonal mean kinetic energy: 



2 



V dx 



and the eddy energy: 



E' = - 
2 



|V^|^ 



dx 



From the zonal mean equation ([5]) we obtain: 



dE 
'dt 



where 



2rrr,E 



V u'Q'dx 



(8) 



(9) 



(10) 



(11) 



is the time rate of change of the zonal mean en- 
ergy due to the eddy induced mean zonal acceleration 
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—u'C'- Similarly, we obtain from the perturbation 
equations (|6al6bp : 



dE' 
~dt' 



-re + r„-r„-r, + F , (12) 



where 



= K u'n'dx , 



(13) 



is the rate of perturbation energy gain due to pertur- 
bation density flux down the mean density gradient. 
This term provides the internal energy source for the 
turbulence. The term 



Ta = a 



n)^dx , 



(14) 



corresponding to resistive coupling is always dissipa- 
tive as is the diffusion: 



The term, F, is the rate of energy input by external 
excitation. This external energy input rate is con- 
stant if the excitation is delta correlated and state 
independent. 

b. The SSST system governing DW-ZF dynamics 

Wc parameterize the nonlinear scattering term ([7]) 
in the eddy equations (pal [6b)) by stochastic forc- 
ing, which is the STM closure (jFarrell and loannou 
1993cl Il994al ). The STM accurately simulates both 
the structure of the eddy field and of the quadratic 
fluxes in shear turbulence including that of the 
Earth's atmosphere, which is a particularly w ell ob- 
served turbulent medium ( FarrcU and loannou 1995 



DelSole 1996 



19981 



Whitaker and Sardeshmukh 
Zhang and Heidi 1 19991 : lDelSolel l2004l) . 

We represent the perturbation fields using Fourier 
components in the poloidal direction, y: 



^yjkix,t)e'''y , n' = ^nfe(x,t)e'^^ , (16) 

fc 



fc 



and discretize the equations in the radial direction, x, 
so that the state tpk = [<l>k,nk]'^ is prescribed by the 
values, for each Fourier component, of the electro- 
static potential and the perturbation density on an 
equally spaced grid. Under the simplifying assump- 
tion that the stochastic forcing has sufficiently short 
temporal correlation that it can be approximated as 
white noise, the second moment statistics of the fiuc- 
tuating field, ■0^, are fully described by the covariance 
matrix Ck =< tpk'4'k > ( <'> denotes ensemble av- 
eraging) which evolves according to the deterministic 



Lyapunov equation: 
dCk 



dt 



Akiv)Ck 



CkKiv) 



Qfc , (17) 



= / " (|VCf + |Vnf)dx . (15) [ 

"'0 



in which Qk is the covariance representing the ensem- 
ble average distribut ion of the stochastic forci ng in 
the radial direction ( Farrell and loannoul 1996h and 
Ak{v) is the linear operator in ([6al I6b|l which de- 
pends affinely on the zonal flow v[x,t). If Q/j repre- 
sents scattering by the advective nonlinearity rather 
than external sources of excitation then a dissipation 
can be added to the linear operators to ensure that 
no net energy is introduced into the system (because 
the nonlinear terms only redistribute energy). Also, 
Qfe can be made an appropriate function of the am- 
plitude of the perturbation variance in order to accu- 
rately parameterize the quadratic nonlinearity of the 
advective Jacobian. More comprehensive cl osures of 
this sort have been used in other contexts (iDelSolel 
2001bHFarren and Ioannoull2009bl) : however, it is suf- 
ficient for our present purposes to use the simplest pa- 
rameterization in which the system is stochastically 
excited with state independent forcing and the be- 
havior of the system is investigated as a function of 
the amplitude of this excitation. 

The Lyapunov equation (|17p determines Ck and 
this covariance in turn determines the ensemble mean 
vorticity flux: 



1, 



ik ( (hk^k 



k^k(Pk 



= 5]|3[diag(C,At) 



(18) 



However, it is the zonal mean vorticity flux that 
appears in the zonal flow equation ([5]) but under the 
ergodic assumption the zonal mean can be replaced 
by the ensemble mean: 



(u'Acf)') = u'A<j>' 



(19) 



This requires that there be many independent real- 
izations of eddy activity in the poloidal direction and 
in that limit we obtain the ensemble mean equations: 

dtv = -Y, [diag(CfeAi)] - r„iJ (20a) 
fc 

^ = Ak{v)Ck +CkAl(v) + Qk . (20b) 



The equation for the turbulence covariance, (|20b[) . 
and the equation for the mean zonal flow, (|20a[) . to- 
gether comprise a closed system for the evolution of 
the zonal flow under the influence of its consistent 
field of turbulent eddies. Although the effects of the 
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ensemble mean turbulent fluxes are retained in this 
system, the fluctuations associated with the turbu- 
lent eddy dynamics are suppressed and the dynamics 
becomes autonomous and deterministic. These SSST 
equations can be interpreted as the dynamical equa- 
tions for the evolution of the quadratic (Gaussian) ap- 
proximation to the ensemble mean probability distri- 
bution of the turbulent DW-ZF system. This concept 
invites novel perspectives such as that of chaos of the 
ensemble mean state of a turbulent system as distinct 
from chaos of a realization of the system. We show 
examples of ensemble mean chaos in DW-ZF turbu- 
lence below. However, the SSST system trajectory is 
often not chaotic but instead asymptotes to a flxed 
point equilibrium and in these cases the dynamics of 
DW-ZF equilibria emerge with great clarity in the 
SSST system. As another illustration of the insight 
provided by this system we note that zonal jets arise 
in SSST as easily analyzed linear instabilities. This 
jet forming instability is an example of a new class 
of instability in fluid dynamics; it is an emergent in- 
stability that arises essentially from the interaction 
between the zonal flow and the turbulence. 

c. Parameters 

Unless otherwise indicated calculations were per- 
formed with 64 points in the x direction and 8 har- 
monics in the (y) direction comprising wavenumbers 
k = [A;o, 3A;o, 5fco, 7A;o, 9fco, llfco, 13fco, 15A:o] with ko = 
0.15 in a doubly periodic channel with Ly = 27r/fco 
and Lx = Ly/A. The stochastic forcing is taken to 
have an identity covariance in vorticity correspond- 
ing to a one grid point correlation, and is normalized 
so that the energy input by the stochastic forcing is 
the same for all zonal wavenumbers. The excitation 
of the electrostatic fleld and the density fleld is corre- 
lated in order to facilitate the adjustment of the two 
fields (similar results are obtained using uncorrelated 
forcing). The amplitude of the stochastic forcing is 
given in terms of the equivalent Urms velocity that 
would be maintained by the forcing with no zonal 
flow and with k = 0. Dissipation parameters used 
are p = IQ-^, a = 1 and < r„ < IQ-*. 

3. Results 

a. Formation of zonal jets starting from a non- 
equilibrium state 

The starting point for a systematic investigation 
of DW-ZF dynamics is the nonlinear SSST system 
initiated in a state lying on its attractor. However, 
the system is commonly thought of as being initiated 
far from its attractor in a state of high turbulence 



intensity but without the corresponding finite am- 
plitude equilibrium jet. There then ensues a rapid 
adjustment process in which the system builds a jet 
corresponding to the turbulence and in the process 
places the system on the SSST attractor. In order 
to study this adjustment process consider the exam- 
ple of the turbulence field associated with a single 
poloidal wavenumber in equilibrium with a strong 
stochastic excitation but without its consistent zonal 
jet. The turbulence field is that obtained from the 
stochastically excited STM but without coupling to 
the zonal flow equation. If the zonal flow equation 
is coupled to the STM at this point to form the in- 
teractive SSST system there ensues rapid formation 
of a consistent zonal jet. We show this rapid devel- 
opment of a jet from a strong initial turbulent field 
at the single global wavenumber to = 5 in Fig. [1] for 
K = and no stochastic excitation. We also show 
the unstable case with k = 1 and in both cases the 
development of the zonal jet is rapid because of the 
feedback between the eddies and the growing zonal 
jet. If as an experiment the eddies are required to 
develop on a fixed jet structure that is not continu- 
ously modified by their dynamics then the resulting 
fluxes build the jet much more slowly revealing that 
rapid jet formation is due to the cooperative DW-ZF 
interaction. The build up of the jet and the sub- 
sequent suppression of the eddy energy occurs due 
to shearing of the eddy field by t he jet, a process 



discussed bv lDiamond et ah ( 20051) and that is seen 



both in simulations (jNumata et aL 2007 ) and obser- 



vations. Because the vorticity flux is proportional to 
the shear ( Farrell and Ioannoulll993b . 2009a) the rate 
of increase of the shear is proportional to the shear 
and the eddy variance, so that if the eddy variance is 
constant exponential growth of the shear results but 
if the eddy variance is also growing during this phase 
an exponential growth ensues with time increasing 
exponent. 

Similar development occurs when there is stochas- 
tic forcing and as a result the turbulence has a full 
spectrum. An example with k = 1 of jet emergence 
from small amplitude random initial conditions in an 
unstable flow with substantial stochastic excitation is 
shown in Fig. [2] and the process of its approach to 
equilibrium is shown in Fig. [3l The eddy induced 
zonal acceleration reaches its peak during this initial 
development (cf. Fig. [2ji). The rapid suppression 
of the eddy variance (cf. Fig. [2t) is caused by en- 
ergy transfer to the zonal flow and by increased dis- 
sipation, Fq, due to increased disequilibrium of the 
electrostatic fleld </>' and the perturbation ion density 
fluctuations n' [cf. the energetics equation p2)) ]. 

After the initial development of the zonal jet by the 
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mechanism of anti-diffusive shear momentum trans- 
port as described above there follows a period of ad- 
justment in which the SSST system attempts to stabi- 
lize the zonal flow and to establish, if the parameters 
allow it, a steady state equilibrium corresponding to a 
fixed point of the SSST equations. This stabilization 
process is shown in Fig. [3] as it sequentially stabilizes 
the perturbation operator A^. As the jet adjusts to 
equilibrium during this phase the flow is dominated 
by large structures and the adjustment has a full wave 
modal character unlike during the initial period of jet 
formation from a state far removed from the system 
attractor in which the dynamics is shear wave anti- 
diffusion dominated being associated essentially with 
rapid distortion of the initial perturbation field. 

b. Structural instability of the zero zonal flow state 
as a function of the amplitude of the stochastic 
excitation in the absence of drifl wave instability, 

K = 0. 

We turn now to dynamics on the attractor of the 
SSST DW-ZF system and first study the case k = 
in which there is no drift wave instability and eddy 
variance is maintained solely by external excitation. 
In the absence of zonal flow the SSST equations (|20al 
I20bp are translationally invariant in the radial direc- 
tion and the vorticity flux m'C vanishes and as a re- 
sult the zero zonal flow is an equilibrium of the SSST 
equations for any stochastic excitation and associ- 
ated turbulence level. The SSST equations can be 
linearized about this zero state ve = and the eddy 
covariance that corresponds to a chosen stochastic 
excitation of this zero state, C^^;, obtained from the 
steady state Lyapunov equation. About this state 
perturbation equations can be obtained for the per- 
turbation zonal velocity, 5v, and perturbation eddy 
covariances, SCk, in the form: 



Sv 
SCk 



L{VE, CkE) 



Sv 
SCk 



(21) 



The growth rate and structure of the most rapidly 
growing eigenmode of L provides insight into the 
mechanism of zonal jet emergence and equil ibration 
in turbulence! Farrell and Ioannoul2003 . 2007 ). Zonal 
jets arise as finite amplitude nonlinear equilibria pro- 
ceeding from the most rapidly growing eigenmode of 
L linearized about the zero state. Note that this jet 
forming instability docs not in general coincide with 
loss of stability of the A^ operators which determine 
the stability of a finite amplitude zonal flow to eddy 
perturbation. 

The SSST system can be linearized about finite 
amplitude SSST equilibria as well as about the zero 



state and the bifurcation structure about these fi- 
nite equilibria can be examined as a function of pa- 
rameters to determine e.g. the circumstances under 
which jet breakdown occurs. It should be noted in 
this context that equilibria of the SSST system are 
necessarily perturbation stable. Consider as an ex- 
ample the SSST stability of the zero zonal flow state, 
[CkEjVE = 0], with K = 0. In this case the zonal 
jet emerges as increase in stochastic excitation, Qfc, 
causes the turbulence level to exceed a threshold at 
which point L becomes SSST unstable. As the am- 
plitude of the excitation, Q^, is increased further this 
bifurcation connects to finite amplitude equilibria in 
which the eddies maintain finite amplitude zonal jets. 
The bifurcation diagram of this zonal flow as a func- 
tion of excitation amplitude is shown in Fig. [4K,c to- 
gether with the associated non-linearly equilibrated 
zonal jets. For weakly supercritical excitation the 
structure of the zonal flow is nearly that of the most 
unstable mode of the L operator but as the excitation 
increases the velocity of the zonal flow asymptotes to 
a constant structure as shown in Fig. [Hi. 

We can understand an important aspect of the 
dynamics of this asymptotic structure by noting 
that as the stochastic excitation increases the zonal 
flow acceleration associated with the ensemble mean 
Reynolds stress divergence: 



< u'C >-- 



UkCl 



(22) 



is comprised of a sum of low zonal wavenumber fluxes 
that decelerate the jet and high wavenumber fluxes 
that accelerate the jet. As excitation and turbu- 
lence level increase the vorticity flux of each com- 
ponent of this sum increases while the sum tends to 
the small residual required to balance the zonal flow 
dissipation because the low wavenumber downgradi- 
ent and high wave number upgradient contrib utions 
very nearly cancel (jFarrell and Ioannoull2009al ) . This 
dynamic can be seen in Fig. [5] in which the struc- 
ture of the vorticity fluxes associated with the equi- 
librium jet in Fig. HJl is shown. In Fig. [5^ it is seen 
that wavenumbers m — 1,3,5 oppose the jet and 
nearly cancel the upgradient contribution from the 
higher wavenumbers. This cancellation becomes all 
the more complete as the excitation increases and the 
equilibrium zonal flow assumes asymptotic form. Be- 
cause the total vorticity flux vanishes in the collision- 
less limit, r„i ~ 0, these equilibria are also the equi- 
libria in this limit (as shown in Fig. ^jp). This demon- 
strates that in turbulence with vanishing coUisional 
damping of the zonal flow there are non-vanishing 
equilibria that are independent of the turbulence in- 
tensity and have the universal structure shown in Fig. 
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I3J1. It should be noted that while this asymptotic 
zonal flow does not depend on the turbulence inten- 
sity for a fixed spectrum of excitation it is sensitive to 
the spectral distribution because the fluxes are upgra- 
dient for high wavenumbers and downgradient for low 
wavenumbers. For example, if only low wavenumbers 
are excited no finite equilibria result for any r,„ as 
all fluxes oppose the jet. Conversely, if only the high 
wavenumbers are excited equilibria arise for r„ 7^ 
associated with upgradient fluxes but in this case the 
equilibrium zonal flow increases secularly with exci- 
tation increase until the jet became structurally un- 
stable. 

We note, in the examples shown and in agreement 
with observations and simulations, that the kinetic 
energy is concentrated in the energy of the zonal jet 
while the eddy kinetic energy is greatly suppressed 
(cf. Fig. ED. 

c. Zonal flow equilibria as a function of the ampli- 
tude of stochastic excitation in the presence of 
drift wave instability, k > 0. 

Introduction of unstable density stratification k > 
makes the zero state perturbation unstable and nec- 
essarily structurally unstable for any stochastic exci- 
tation. These unstable eddy perturbations augment 
the turbulence and facilitate formation of zonal jets. 
An example with k = 1 of jet emergence from small 
amplitude random initial conditions in an unstable 
flow with substantial stochastic excitation are shown 
in Fig. m and Fig. H 

The radial distribution at various poloidal 
wavenumbers of the equilibrium particle flux and of 
the acceleration by the Reynolds stress are shown in 
Fig. [7^,c. The eddy induced acceleration of the zonal 
flow by small wavenumber eddies is downgradient, 
opposing the jet, while the acceleration due to larger 
wavenumbers is upgradient, as for the case k = 0. 
This cancellation implies, as for the case with k — 0, 
that the equilibrium flow asymptotes to a fixed struc- 
ture as the amplitude of the forcing increases and that 
this asymptotic flow is also the equilibrium flow in 
the collisionless limit, = 0. Similar equilibria with 
zero coUisional dar nping have also b een seen in turbu- 



lence simulations ( Lin et al. 20001 ). The asymptotic 



equilibrium flow shown in Fig. [5] is found to depend 
only weakly on k. For = this is the universal 
equilibrium flow for all forcing amplitudes and for all 
K. However this equilibrium is structurally unstable 
for large values of will be discussed. 

The eddy kinetic energy peaks at the gravest 
poloidal scale, m — I. It is important to note that 
it is at large scales that most of the eddy energy re- 



sides and also it is the large scales that are respon- 
sible for the particle flux (the particle flux peaks for 
TO = 3 as shown in Fig. Uh) as is a lso found in tur- 
bulent simulations ( Bos et aDl2008h . The dominance 
of large scales in the eddy variance and fluxes is con- 
sistent with these scales being the least damped (cf 
Fig. [3]), however the eddy structure does not assume 
the structure of the least damped modes. We show 
the structure of the least damped mode and the dis- 
tinct structure of the top EOF of the covariance ma- 
trix (this is the eigenfunction associated to the largest 
eigenvalue of C^) for the gravest poloidal scale to = 1 
in Fig. [9^,f. We also show the first stochastic optimal 
which is the structure of the excitation that would 
produce, if the flow were forced stochastically with 
this structure at unit amplitud e, the highest eddy 
energ y at statistical equilibrium (^ Farrell and loannoul 
19961 ). The difference in the structure of the top EOF, 



the least damped mode and the stochastic optimal re- 
veal the degree of non-orthogonality of the modes of 
the operator which is related to their non-normality 
as these would be identical if the system we re normal 
( Farrell and Ioannoulll994b : loannou 1995 ). 

The non-normality of the H-W system is central to 
its dynamics. In order to appreciate its role consider 
the frequency spectrum of the total eddy variance re- 
sulting from excitation unbiased in time and space of 
the linearized H-W equations. This can be obtained 
by Fourier transforming the perturbation equations 
al I6b[) written in the form: 



dt 



Aki'k + Fi] 



(23) 



where 77 is Gaussian white noise and F gives the radial 
structure of the forcing related to the noise covariance 
Qk in ini) by Qfe = FF^, to obtain: 



Muj) = R{oj)Ffi{uj) , 



(24) 



where variables that depend on lo denote the Fourier 
amplitudes, i.e., 

M^) = TT / ^fcWe-'"' dt , (25) 

and fi{uj) is the Fourier amplitude of the Gaussian 
noise. The resolvent R(a;) determines the structure 
of the response and is given by 

Rfc(tj) = (ic^I- Afe)-i . (26) 

We form the correlation matrix 

Ck{uj) = (m^)M^A = Rfe(w)QfcRfe(w)t (27) 
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and proceed to calculate the perturbation energy 
power spectrum as a function of phase velocity as 



tion by small eddies as would be the case if it were 
diffusive. 



E{c) 



k 



trace 



1/2 



(28) 



and Mfc is the energy metric defined so that Ek = 
iJ^llsAkipk is the perturbation energy. The power spec- 
trum is shown in Fig. [10] both for /t = and n — 1 
along with the equivalent normal response which is 
obtained by calculating the power spectrum by re- 
placing Afc by a diagonal matrix with elements its 
eigenvalues. If the forcing covariance were the iden- 
tity the equivalent normal response would be given 
by the resonance formula: 



E 



1 



lU! ■ 



in 



kjl 



(29) 



where iflk are the eigenvalues of Afc. The equivalent 
normal power spectrum is equal to the power spec- 
trum when Afc is a normal matrix and the forcing 
is an identity, otherwise the power spectrum exceeds 
the equivalent normal power spectrum and the dif- 
ference reveals the degree of non-normality. The dif- 
ference reflects the excess power that is maintained 
by the system against friction because of the non 



orthogonality of the eigenmodes (jFarrell and loannou 
[l994b; l oannou J995i ). Note that the power peaks 
at phase speeds near the maximum and minimum 
velocity of the zonal flow. An asymmetry develops 
as K increases with power becoming concentrated in 
the prograde jet reflecting the increased instability of 
the prograde jet as compared with the retrograde jet 
when K > 0. Because the frequency response arises 
primarily from the gravest poloidal wavenumber this 
double peak in the turbulence spectrum as a func- 
tion of phase speed is reflected in the frequency spec- 
trum with a double peak at w = fc- 



mtn ^max i 



where 

kmin is the poloidal wavenumber corresponding to 
the gravest mode and Vmax is the maximum velocity 
of the zonal flow. Similar strongly peaked spectra in- 
dicative of coherent large scale structures in zon al jet 
equilibria have been observed (iBush et al.ll200,'^ . 

The particle flux at equilibrium reflects the struc- 
tures producing it. This flux reaches a maximum as 
a function of poloidal wavenumber at m = 3 as seen 
in Fig. [7|d. The flux is downgradient where the jet is 
prograde and becomes upgradient where the jet is ret- 
rograde. The difference between the upgradient and 
downgradient particle fluxes leads to a small down- 
gradient residual which is responsible for the eddy 
energy source. The regions of upgradient flux show 
that the particle flux is produced by large coherent 
structures rather than resulting from random advec- 



d. Zonal flow equilibria for k > 

The dependence of zonal flow equilibria on the 
amplitude of the stochastic excitation in the pres- 
ence of an internal energy source (k = 1) is sim- 
ilar to that of zonal flows in the case without an 
internal energy source (k = 0). We find equilibria 
in the coUisionless limit, r„j — 0, and these exist 
for all forcing amplitudes. These equilibria are in- 
dicated by a dashed line in the bifurcation diagram 
in Fig. [TTt along with the equilibria that result for 
r,„ = 10""*. The equilibria for non zero damping tend 
to the equilibria for = as the stochastic excita- 
tion increases. This asymptotic is reflected in the 
eddy induced zonal flow acceleration which asymp- 
totes as the stochastic excitation increases (shown in 
Fig. I lib ). The eddy kinetic energy at equilibrium 
increases with the amplitude of the stochastic exci- 
tation and is minimized for zero coUisional damping, 
r,„ = 0. The particle flux, measured by the average 
value Tn/Lx, increases with stochastic excitation and 
for zero mean coUisional damping the particle flux 
is increasing quadratically with stochastic excitation 
according to Tn/Lx = 0.0265u^„s. From this it is 
clear that it would be desirable to operate a device 
at low stochastic excitation levels and with reduced 
mean coUisional damping if maximizing confinement 
is the goal. All the equilibria of Fig. [Tl] are struc- 
turally stable for the chosen parameters. However the 
basin of attraction of the equilibria is not the whole 
space. Also note that there are no equilibria with 
rm = 10~^ for stochastic excitations smaller than 
Urms = 0.065. 

Stochastic excitation, which augments the turbu- 
lence, is important for the equilibration process. In 
the absence of stochastic excitation the eddy field 
is dominated by the fastest growing modes and the 
structure of the covariance is not of high enough rank 
to comprise the diversity of structures required to 
produce equilibration. At zero or very low stochastic 
excitation a vacillation regime is found as occurs for 
slightly supercritical states in baroclinic turbulence 
(jPed loskv 1977) while for sufficiently high excitation 
and associated turbulence levels one obtains equilib- 
ria. These equilibria for substantial stochastic excita- 
tion (i.e. Urms > 0.1) are not only structurally stable 
but also have a basin of attraction that spans the 
whole space. However, as the excitation and the sup- 
ported turbulence is reduced the basin of attraction 
of the equilibria shrinks and finally at a critical value 
equilibria cease to exist. Operationally, states with 
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low stochastic excitation and small particle fluxes can 
be approached by first obtaining an equilibrium by in- 
creasing the stochastic excitation and then adiabati- 
cally adjusting the parameters to reach these isolated 
in parameter space states. 

The vacillation regime mentioned above is not a 
vacillation of the trajectory of a single realization of 
the turbulence but rather a vacillation regime in the 
trajectory of the probability density function of the 
turbulence in the (Gaussian) SSST approximation. 

We show in Fig. [THa state of chaotic DW-ZF fluc- 
tuations at K = 1 with very weak forcing [producing 
equivalent u^ms — O(10~^) /{ps^d) ] a-nd zero mean 
collisional damping. We have determined that there 
exists an equilibrium state but this equilibrium state 
has a small basin of attraction and can not be ap- 
proached from the SSST initial conditions chosen in 
this example (which are low turbulence levels, and 
very small zonal flow). A similar chaotic state per- 
sists for these parameters when the collisional damp- 
ing is raised to = 10~^, but for that damping 
there is no SSST equilibrium underlying this state 
(cf Fig. [TTbV In Fig. WBn , we see the initial devel- 
opment of the zonal flow, followed by an adjustment 
period, but unlike the case with strong stochastic ex- 
citation shown in Fig. [2l which adjusted to equi- 
librium by stabilizing the perturbations, the insta- 
bility remains and alternating periods ensue of high 
eddy activity (low zonal flow) and low eddy activity 
(high zonal flow). The fluctuations settle to a chaotic 
bursting pattern in the zonal flow and the eddy vari- 
ables as shown in Fig. [131 The eddy variables, the 
particle flux at a specific location, the eddy kinetic 
energy and the integrated particle flux, have a saw- 
tooth structure in which a slow build up of the eddy 
variance associated with the underlying instability is 
followed by a rapid collapse of the eddy flelds as the 
zonal flow develops and converts the eddy energy to 
mean zonal energy over an advective time scale. The 
mean zonal kinetic energy exhibits a sawtooth be- 
havior in which the mean develops very rapidly and 
then slowly adjusts under the influence of the weak in- 
duced mean eddy accelerations. Such sawtooth struc- 
tures have been commonly observe d and simulated 
(|Wagneill2007l : iDiamond et al.ll2005h . 

For the same parameters that we obtained the 
chaotic regimes shown in Fig. [13] there exists an iso- 
lated stable equilibrium with a limited basin of at- 
traction. This equilibrium state can be elicited by 
impulsively introducing any SSST zonal flow that is 
stable at these parameters. Immediately upon intro- 
duction of the zonal flow the eddy energy and the 
particle flux are quenched and the flow asymptoti- 
cally relax to the equilibrium flow as shown in Fig. 



1141 If the parameters do not support an equilibrium 
DW-ZF state then the zonal flow eventually breaks 
down and a chaotic regime ensues. For example if the 



collisional damping is raised to r,„ = 10 



there 



is no equilibrium at this amplitude of stochastic ex- 
citation and in time 0(l/7',„) the imposed zonal flow 
reverts to a chaotic state. 

Regime transition can be controlled using stochas- 
tic excitation as a control parameter. As the exci- 
tation increases the chaotic bursting gives way to a 
quasi-periodic regime and by further increasing the 
stochastic excitation a fixed point DW-ZF equilib- 
rium jet state is established as shown in Fig. [T5l 
Having obtained an equilibrium jet state we then re- 
duce the stochastic excitation (shown in Fig. fTH]) 
and find that the jet persists as the stochastic ex- 
citation is reduced and both the eddy kinetic energy 
and the particle fiux vanish with the excitation. This 
equilibrium state exists at the same parameter val- 
ues for which periodic and chaotic behavior are ob- 
tained. Hysteretic transition between a steady zonal 
flow state and a chaotic turbulent state is common 
in turbulent systems such as sheared boundary layer 
flows which exist in laminar and turbulent states at 
the same parameter values. 

Dependence of zonal flow, eddy variances and 
fluxes at equilibrium on mean collisional damping is 
shown Fig. 1171 the particle flux increases with mean 
collisional damping, as does the eddy energy while 
the zonal flow velocity decreases, as is al so found in 
turbulence simulations ( Itoh et al. 2007bf l. 



e. Loss of structural stability at large k 

We now investigate zonal flow equilibria as a func- 
tion of K. These equilibria, as already discussed, 
are most easily initialized at high stochastic exci- 
tation amplitude and low mean collisional damp- 
ing. We study the dependence of these equihb- 
ria on k at high turbulence levels (with equivalent 
Urms — 0.34 /{psUci))- The maximum zonal flow 
speed is shown in Fig. \18h as a function of k and 
the mean particle flux averaged over the whole do- 
main is shown in Fig. 118b . The particle flux is seen 
to initially increase linearly with n. The equilibria 
are globally attracting up to about k — 1.5, for the 
parameters of this problem, but the basin of attrac- 
tion contracts as k is increased until the flow becomes 
structurally unstable at k = 2.534, and no equilibria 
exist for larger values of k. Although equilibria exist 
for K > 1.5 these equilibria can not be reached from 
the above listed flxed parameters starting from any 
initial condition, but they can be reached by flrst es- 
tablishing an equilibrated state at a lower value of 
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K and then increasing k adiabatically; although op- 
erationally these states are most readily established 
by first going to higher stochastic excitation, corre- 
sponding to a higher level of turbulence, then increas- 
ing K and finally reducing the excitation. 

The equilibrated flows and the corresponding min- 
imum damping decay rate of the least damped mode 
at each poloidal wavenumber, m are shown in Fig. 
[in] for the critical Kc — 2.534 and for the smaller un- 
stable stratifications n = 2.52 and k = 2.0425. As 
the critical value of is approached the poloidal, 
m — 5 wave, tends towards instability. However, 
as — ^ the damping decay rate of the least 
damped mode approaches, for the chosen parameters 
kcimax — 0.12a;ci whilc the fluxes and the equilib- 
rium zonal flow tend to diverge as Kc is approached 
and no equilibrium flows can be sustained for k > Kc 
and transition to a time varying state occurs. This 
result shows that the jet first loses structural stability 
as a function of k rather than modal stability. 

4. Discussion 

There are a number of points we wish to emphasize 
in connection with the above results: 

1. A novel concept arising from SSST is that of 
the structural stability boundary for zonal flow 
breakdown as distinct from breakdown related 
to shear instability of the zonal flow. 

2. Multiple DW-ZF regimes are predicted to ex- 
ist in parameter space including a regime of 
steady zonal flows as well as regimes of pe- 
riodic, quasi-periodic and chaotic bursting or 
"sawtoothbehavior. These regimes provide op- 
portunity for placing and manipulating conflne- 
ment devices to be in a desired dynamical state 
between high and low conflnement. 

3. SSST predicts that isolated DW-ZF equihbria at 
high K are not connected continuously to lower k 
states but that these states can be reached either 
using external turbulence excitation or flnite am- 
plitude state perturbation to promote the system 
between these equilibria. 

4. A mechanism for introducing and modulating 
turbulence levels is predicted to provide a pow- 
erful control parameter for placing the DW-ZF 
system in desired confinement states. 

5. In the limit of vanishing zonal flow coUisional 
damping a universal DW-ZF state is supported 
in which a precise balance between down gra- 
dient momentum transport by small wavenum- 
bers and upgradient transport by high zonal 



wavenumbers occurs. This asymptotic equilib- 
rium predicts that band limiting turbulence can 
prevent formation of stable equilibrium zonal 
flows. 

6. Density fluxes are not diffusive but rather are 
primarily produced by large scale structures. 
Robust fluxes both up and down the mean gra- 
dient occur and it follows that particle transport 
analysis requires a full wave solution. 

5. Conclusion 

Emergence of zonal jets in turbulent flow and the 
relation of these jets to the statistical equilibrium of 
the turbulent state is a problem of great theoreti- 
cal and practical interest. This problem is partic- 
ularly compelling in the case of turbulent plasmas 
because of the relationship of zonal jets to the H 
states that limit turbulent transport of particles and 
heat in magnetic confinement fusion devices. DW- 
ZF interaction dynamics is responsible for the gen- 
eration and regulation of these zonal flows so it fol- 
lows that prospects for predicting and controlling the 
H state require improvement in fundamental under- 
standing of the mechanism underlying the statistical 
steady state of zonal jets in drift wave turbulence. 
In this work we applied the methods of SSST to the 
Hasegawa-Wakatani model to study the emergence, 
stability and effect on transport of zonal jets in the 
DW-ZF system. We find robust zonal jet formation 
in agreement with both experiment and simulation 
and obtain parameter requirements for jet formation 
and breakdown. We find multiple regimes including 
chaotic, periodic and steady and show that externally 
imposed turbulence and finite amplitude zonal flow 
perturbations can be used to control regime transi- 
tion. We find suppression of particle transport by 
zonal fiows and show that this transport is not diffu- 
sive in nature. These results provide a basis for pre- 
diction and controlling confinement regimes in DW- 
ZF turbulence. 
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1 Initial jet formation by the rapid adjustment process starting from a state of strong turbulence 
for the cases (a) k = (no instability) and (b) n — 1 (strong instability). Shown are eddy 
kinetic energy (dashed) and mean zonal kinetic energy (solid) as a function of time. The eddy 
field is limited to global zonal wavenumber m = 5 and there is no stochastic excitation. . . 

2 Transient development of an equilibrium zonal jet. (a) time development of the mean ki- 
netic energy of the zonal flow, Em/ (noTe) (solid), the mean eddy kinetic energy Ke/{noTe) 
(dashed) and the total particle flux over the channel TnUJci/inoTe) (dash-dot), (b): zonal 
velocity, V/{psUJci), as a function of the radial direction and time, (c): eddy kinetic energy, 
\ogiQ{Ke/ (npTe)), as a function of the radial direction and time, (d): eddy induced zonal flow 
acceleration, — < u'C' > /(psW^J, as a function of the radial direction and time. The param- 
eters are: k — 1, r„i — 10~'^Uci and the stochastic excitation has equivalent r.m.s. velocity of 

O.MpsLOci 

3 Evolution of the zonal flow and its associated spectrum for the example in Fig. 2. Left panels: 
zonal flow structure at T = 10,30 and at equilibrium. Center panels: spectrum {cr,kci) of 
Afe for the flow in the corresponding panel for zonal wavenumber to = 3. The continuous line 
indicates the velocity interval spanned by the zonal flow. At equilibrium the instabilities have 
been stabilized. Right panels: the largest growth rate for a given zonal wavenumber, to. At 
equilibrium the least stable mode corresponds to the gravest zonal wavenumber 

4 (a) Maximum zonal flow velocity as a function of stochastic excitation for k — 0. Stochastic 
excitation is measured by the Urms that would have been maintained in the absence flow. 
For the chosen parameters the critical forcing required to form zonal flows is Urms = 7.8 x 
10^^ /{psUJci)- (b) Corresponding equilibrium zonal zonal flows: the larger velocity corresponds 
to forcing denoted with a circle in panel (a), while the smaller velocity corresponds to the 
parameters denoted with a square in (a), (c) Continuation of the bifurcation diagram of (a) 
to larger forcing values. Note that as the forcing increases the maximum zonal flow velocity 
asymptotes to a constant, (d) The asymptotic zonal flow at large forcing. The coUisional 
damping of the mean is — 10~^(jJci 

5 Structure of the eddy induced zonal flow acceleration — < u'C > as a function of radius. The 
solid line is the total flux summed over all zonal wavenumbers multiplied by 100 (at equilibrium 
this is equal to 100r,„w). The dashed line is the acceleration induced by wavenumbers to = 
7 — 15. These higher waves build the zonal flow. The dash-dot line is the acceleration induced 
by the small wavenumbers to = 1 — 5 which tend to destroy the zonal flow. Left panel: for 
mean collisional damping rm — 10^** and the equilibrium flow in Fig. [TJi. Right panel: For 
r„i = (here the cancellation between downgradient and upgradient fluxes is perfect). . . . 

6 Ratio of mean zonal kinetic energy to eddy kinetic energy as a function of stochastic excitation 
(solid). Ratio of the mean zonal kinetic energy to the eddy kinetic that would have been 
maintained in the absence of the zonal flow (dashed). For small excitations there is no zonal 
flow and the ratio vanishes, also for large excitations the flow asymptotes to a constant and 
again the ratio vanishes. For intermediate excitations the zonal flow energy is two to three 
orders of magnitude larger and the turbulence energy is dominated by the zonal flow energy. 
The zonal flow suppresses the eddy energy by approximately an order of magnitude. For k = 
and r„j — lO^^cJci 
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7 (a) Structure in radius of the particle flux at equilibrium. The particle flux is not diffusive, as 
it has a distinct structure and there is a region of upgradient flux that would correspond to 
a negative diffusion coefficient, (b) The integrated particle flux at equilibrium as a function 
of zonal wavenumber, to. (c) The structure of the eddy acceleration — < u' C,' > produced 
by the zonal modes. The thick solid line is the total vorticity flux which maintains the zonal 
flow against dissipation shown in Fig. [51 The opposing fluxes (solid and dashed) is the 
flux associated with wavenumbers to = 1, 3 while the supporting fluxes (solid and dashed-dot) 
correspond to the higher wavenumbers m = 5,7. (d) The energy of the eddy fleld as a function 
of zonal wavenumber. The eddy kinetic energy peaks at the gravest zonal mode to = 1. The 
case is for k = 1 = 10^^ ujd and stochastic excitation equivalent to r.m.s. velocity of 
O.SApsUJci [23 

8 Zonal flow at equilibrium as a function of radius. Dashed: with no coUisional damping of the 
mean (r^ = 0); solid: with r™ — lO^^ujci- The case is for k — 1, and stochastic forcing with 
equivalent r.m.s. velocity of 0.34psWcj \24 

9 (Color online) Top row: the top EOF of the eddy covariance of the component of the eddy 
fleld with zonal wavenumber to 1 (on the left: perturbation electric field, on the right: 
perturbation density) . The first EOF accounts for 32% of the total energy of the eddy field 
at this wavenumber. Middle row: the stochastic optimal. The stochastic optimal produces 
20% of the eddy energy at this wavenumber. Bottom row: the least stable eigenvalue of the 
operator at to = 1. The associated growth rate is ka = —O.WuJci. For the equilibrium zonal 
flow obtained at k — 1 with stochastic excitation equivalent to equivalent r.m.s. velocity of 
O.SApsUci [25 

10 (Color online) Power spectrum of the eddy energy as a function of phase speed Cr (solid). The 
dashed line is the equivalent normal response and circles mark the maximum and minimum 
velocity of the equilibrium flow, (a) for k = 0. (b) for k — 1. The case is for equivalent r.m.s. 
velocity of 0.34psLOci and — 10~^ [H 

11 (Color online) (a) Particle flux as a function of stochastic excitation measured by equivalent 
Urms', for r„j = lO^^ujci (solid) and for = (dashed), (b) Maximum vorticity flux < uX' > 
as a function of stochastic excitation, (c) Maximum equilibrium zonal flow velocity as a 
function of stochastic excitation; for — lO"'*^^,^ (solid) and for — (dashed), (d) 
Mean eddy kinetic energy as a function of stochastic excitation. Also shown is the eddy 
kinetic energy maintained against dissipation in the absence of flow as a function of stochastic 
excitation (dashed-dot). For k = 1 [27 

12 A chaotic state (analysis of perturbed trajectory differences reveals this system to be chaotic 
with Lyapunov exponent 0.02a;ci)- i^)- Zonal flow energy (solid), and eddy kinetic energy 
(dashed) as a function of time, (b): zonal velocity, V/{psUJci), as a function of radius and 
time, (c): eddy kinetic energy, log]^g(/ire/(noTe)), as a function of radius and time, (d): eddy 
induced zonal flow acceleration, — < u'(^' > /{paUi'^i), as a function of radius and time. The 
parameters are: k = 1, = and the stochastic excitation has equivalent r.m.s. velocity of 
0.34 X 10~'^ psOJci- For these values there exists an equilibrium zonal flow with a limited basin 
of attraction, and this equilibrium state can not be approached from initial states with small 
zonal flows |28 

13 For the case shown in Fig. [T^l (a) particle flux at a single location as a function of time; (b) 
zonal flow kinetic energy; (c) eddy kinetic energy; (d) average particle flux [2^ 

14 A chaotic state is laminarized by impulsive introduction of a stable zonal flow at LUdt = 310. 
The zonal flow subsequently asymptotically approach the equilibrium zonal flow that exists 
for these parameter values, (a): zonal velocity, V/{psLUci), as a function of radius and time, 
(b): Zonal flow energy (solid), and eddy kinetic energy (dashed) as a function of time, (c): 
Mean particle flux as a function of time For k = I and = 

15 A chaotic state (0 < uJdt < 400) becomes quasi periodic (450 < ujdt < 2800) and then settles 
to an equilibrium as stochastic excitation increases (0.34 x lO^^ps^ci < Urms < 0.1823psa;ci). 
(a): zonal velocity, V/ (psUJd), as a function of radius and time, (b): Zonal flow energy (solid), 
and eddy kinetic energy (dashed) as a function of time, (c): Mean particle flux as a function 

of time (d) Stochastic excitation as a function of time. For k = 1 and = |31 
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16 Continuation of Fig. [T5l The stochastic excitation is decreased to its initial value {u„ns = 
0.34 X psUici)- The zonal flow persists while the eddy kinetic energy and the particle flux 
vanish with the excitation, (a): zonal velocity, V/{ps(jJci), as a function of radius and time, 
(b): Zonal flow energy (solid), and eddy kinetic energy (dashed) as a function of time, (c): 
Mean particle flux as a function of time (d) Stochastic excitation as a function of time. For 

K = 1 and Td = |32 

17 Equilibrium state diagnostics as a function of mean coUisional damping, (a) Particle flux, (b) 
Maximum vorticity flux, (c) Maximum equilibrium zonal flow velocity, (d) Mean eddy kinetic 



energy. The case is for k — 1, and stochastic forcing with equivalent r.m.s. velocity of Q.iApgUJci- l33 

18 Equilibrium state diagnostics as a function of density gradient, n (a): Maximum velocity of the 
equilibrium zonal flow, (b) The mean particle flux (solid). The mean particle flux increases 
at first linearly as Q.Qbn/Lx (dashed). The parameters arc: r„i = 10~'' and the stochastic 
excitation supports equivalent r.m.s. velocity of Q.ZApaLOci l34 

19 Approach to structural instability as a function of k. Top: Zonal flow velocities as the critical 
Kc — 2.534 is approached. Bottom: The corresponding maximum growth rate of perturbations 
as a function of poloidal wavenumber, m. Solid: k = 2.534, dash: k — 2.52, dash-dot: for 
K = 2.0425. The parameters are: r,„ — 10"'' and the stochastic excitation has equivalent 
r.m.s. velocity of 0.34ps'^ci l35 
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Fig. 1. Initial jet formation by the rapid adjustment process starting from a state of strong turbulence 
for the cases (a) k = (no instability) and (b) k = 1 (strong instability). Shown are eddy kinetic energy 
(dashed) and mean zonal kinetic energy (solid) as a function of time. The eddy field is limited to global 
zonal wavenumber m = 5 and there is no stochastic excitation. 
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Fig. 2. Transient development of an equilibrium zonal jet. (a) time development of the mean kinetic energy 
of the zonal flow, Em/ (noTe) (solid), the mean eddy kinetic energy K^/ (noTf.) (dashed) and the total particle 
flux over the channel TntOci/inQTe) (dash-dot), (b): zonal velocity, V/ipsCOd), as a function of the radial 
direction and time, (c): eddy kinetic energy, logio(-K'e/('^o7e)), as a function of the radial direction and 
time, (d): eddy induced zonal flow acceleration, — < u'(' > /(psW^j), as a function of the radial direction 
and time. The parameters are: k = 1, = lO^'^ujci and the stochastic excitation has equivalent r.m.s. 
velocity of O.iApsOJci. 
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Fig. 3. Evolution of the zonal flow and its associated spectrum for the example in Fig. 2. Left panels: 

zonal flow structure at T = 10,30 and at equilibrium. Center panels: spectrum {cr,kci) of for the flow 
in the corresponding panel for zonal wavenumber m = 3. The continuous line indicates the velocity interval 
spanned by the zonal flow. At equilibrium the instabilities have been stabilized. Right panels: the largest 
growth rate for a given zonal wavenumber, m. At equilibrium the least stable mode corresponds to the 
gravest zonal wavenumber. 
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Fig. 4. (a) Maximum zonal flow velocity as a function of stochastic excitation for k — 0. Stochastic 
excitation is measured by the Urms that would have been maintained in the absence flow. For the chosen 
parameters the critical forcing required to form zonal flows is Urms = 7.8 x 10~'^ /{psUJci)- (b) Corresponding 
equilibrium zonal zonal flows: the larger velocity corresponds to forcing denoted with a circle in panel (a), 
while the smaller velocity corresponds to the parameters denoted with a square in (a), (c) Continuation 
of the bifurcation diagram of (a) to larger forcing values. Note that as the forcing increases the maximum 
zonal flow velocity asymptotes to a constant, (d) The asymptotic zonal flow at large forcing. The coUisional 
damping of the mean is r,„ — 10~^ujci- 
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Fig. 5. Structure of the eddy induced zonal flow acceleration — < u'C,' > as a function of radius. The solid 
line is the total flux summed over all zonal wavenumbers multiplied by 100 (at equilibrium this is equal to 
IQQTrav). The dashed line is the acceleration induced by wavenumbers m — 7 — 15. These higher waves build 
the zonal flow. The dash-dot line is the acceleration induced by the small wavenumbers to = 1 — 5 which 
tend to destroy the zonal flow. Left panel: for mean coUisional damping rm = 10~* and the equilibrium flow 
in Fig. [TJi. Right panel: For r,n = (here the cancellation between downgradient and upgradient fluxes is 
perfect). 
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Fig. 6. Ratio of mean zonal kinetic energy to eddy kinetic energy as a function of stochastic excitation 
(solid). Ratio of the mean zonal kinetic energy to the eddy kinetic that would have been maintained in 
the absence of the zonal flow (dashed) . For small excitations there is no zonal flow and the ratio vanishes, 
also for large excitations the flow asymptotes to a constant and again the ratio vanishes. For intermediate 
excitations the zonal flow energy is two to three orders of magnitude larger and the turbulence energy is 
dominated by the zonal flow energy. The zonal flow suppresses the eddy energy by approximately an order 



of magnitude. For k — and 
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Fig. 7. (a) Structure in radius of the particle flux at equilibrium. The particle flux is not diffusive, as it has 
a distinct structure and there is a region of upgradient flux that would correspond to a negative diffusion 
coefficient, (b) The integrated particle flux at equilibrium as a function of zonal wavenumber, m. (c) The 
structure of the eddy acceleration — < u'C^' > produced by the zonal modes. The thick solid line is the 
total vorticity flux which maintains the zonal flow against dissipation shown in Fig. [H The opposing fluxes 
(solid and dashed) is the flux associated with wavenumbers to = f , 3 while the supporting fluxes (solid and 
dashed-dot) correspond to the higher wavenumbers m — 5,7. (d) The energy of the eddy field as a function 
of zonal wavenumber. The eddy kinetic energy peaks at the gravest zonal mode to = f . The case is for k ^ 1 
r,„ = fO""* LOci and stochastic excitation equivalent to r.m.s. velocity of O.SApgUJci- 
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Fig. 8. Zonal flow at equilibrium as a function of radius. Dashed: with no coUisional damping of the mean 
{fm — 0); solid: with = lO^^Wci- The case is for k = 1, and stochastic forcing with equivalent r.m.s. 
velocity of Q.'i'ipsLJci- 
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Fig. 9. (Color online) Top row: the top EOF of the eddy covariance of the component of the eddy field 
with zonal wavenumber m — I (on the left: perturbation electric field, on the right: perturbation density). 
The first EOF accounts for 32% of the total energy of the eddy field at this wavenumber. Middle row: the 
stochastic optimal. The stochastic optimal produces 20% of the eddy energy at this wavenumber. Bottom 
row: the least stable eigenvalue of the operator at m = 1. The associated growth rate is kci = — O.lSwci- 
For the equilibrium zonal flow obtained at k = 1 with stochastic excitation equivalent to equivalent r.m.s. 
velocity of O.SApsOJci 
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Fig. 10. (Color online) Power spectrum of the eddy energy as a function of phase speed Cr (solid). The 
dashed line is the equivalent normal response and circles mark the maximum and minimum velocity of the 
equilibrium flow, (a) for k = 0. (b) for k = 1. The case is for equivalent r.m.s. velocity of O.SApsUJci and 
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Fig. 11. (Color online) (a) Particle flux as a function of stochastic excitation measured by equivalent Urms] 
for Tm = lO~'^tiJci (solid) and for r,„ = (dashed), (b) Maximum vorticity flux < u'C > as a function of 
stochastic excitation, (c) Maximum equilibrium zonal flow velocity as a function of stochastic excitation; 
for Tm = 10~^ujci (solid) and for r„i — (dashed), (d) Mean eddy kinetic energy as a function of stochastic 
excitation. Also shown is the eddy kinetic energy maintained against dissipation in the absence of flow as a 
function of stochastic excitation (dashed-dot). For k = 1. 
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Fig. 12. A chaotic state (analysis of perturbed trajectory differences reveals this system to be chaotic 
with Lyapunov exponent 0.02ujci)- (a): Zonal flow energy (solid), and eddy kinetic energy (dashed) as 
a function of time, (b): zonal velocity, V/{psUJci), as a function of radius and time, (c): eddy kinetic 
energy, \ogiQ{Ke/{noTe)), as a function of radius and time, (d): eddy induced zonal flow acceleration, 
— < m'C' > /{Ps^ti), as a function of radius and time. The parameters are: k = 1, = and the stochastic 
excitation has equivalent r.m.s. velocity of 0.34 x psUici- For these values there exists an equilibrium 
zonal flow with a limited basin of attraction, and this equilibrium state can not be approached from initial 
states with small zonal flows. 
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Fig. 14. A chaotic state is laminarized by impulsive introduction of a stable zonal flow at LUdt = 310. The 
zonal flow subsequently asymptotically approach the equilibrium zonal flow that exists for these parameter 
values, (a): zonal velocity, V/{psLUci), as a function of radius and time, (b): Zonal flow energy (soUd), and 
eddy kinetic energy (dashed) as a function of time, (c): Mean particle flux as a function of time For k — 1 
and rd = 0. 
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Fig. 15. A chaotic state (0 < ujcit < 400) becomes quasi periodic (450 < uicit < 2800) and then settles to an 
equihbrium as stochastic excitation increases (0.34 x lO^^psWci < Urms < 0.1823/9sa;ci). (a): zonal velocity, 
V/ {psUJci)j as a function of radius and time, (b): Zonal flow energy (solid), and eddy kinetic energy (dashed) 
as a function of time, (c): Mean particle flux as a function of time (d) Stochastic excitation as a function of 
time. For k = 1 and = 0. 
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Fig. 16. Continuation of Fig. [151 The stochastic excitation is decreased to its initial value (urms — 
0.34 X 10^"^ psOJci)- The zonal flow persists while the eddy kinetic energy and the particle flux vanish with the 
excitation, (a): zonal velocity, V/{psLOci), as a function of radius and time, (b): Zonal flow energy (solid), 
and eddy kinetic energy (dashed) as a function of time, (c): Mean particle flux as a function of time (d) 
Stochastic excitation as a function of time. For k, = 1 and rd = 0. 
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Fig. 17. Equilibrium state diagnostics as a function of mean coUisional damping, (a) Particle flux, (b) 
Maximum vorticity flux, (c) Maximum equilibrium zonal flow velocity, (d) Mean eddy kinetic energy. The 
case is for k = 1, and stochastic forcing with equivalent r.m.s. velocity of O.SApgUici- 
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Fig. 18. Equilibrium state diagnostics as a function of density gradient, k (a): Maximum velocity of the 
equilibrium zonal flow, (b) The mean particle flux (solid). The mean particle flux increases at first linearly 
as 0.05k/ Lx (dashed). The parameters are: = 10"'* and the stochastic excitation supports equivalent 
r.m.s. velocity of O.^ApsLOd- 
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Fig. 19. Approach to structural instability as a function of k. Top: Zonal flow velocities as the critical 
Kc — 2.534 is approached. Bottom: The corresponding maximum growth rate of perturbations as a function 
of poloidal wavenumber, m. Solid: k = 2.534, dash: k — 2.52, dash-dot: for k = 2.0425. The parameters 
are: = 10"'^ and the stochastic excitation has equivalent r.m.s. velocity of O.^ipsWd- 
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